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ABSTRACT 



This FepoFt compares three finite element formulations of the 
linearized shallow-water equations which are applied to the geostrophic 
adjustment process. The three corresponding finite difference schemes are 
also included in the study. The development follows Schoenstadt (1980) 
wherein the spatially discretized equations are Fourier transformed in x, 
and then solved with arbitrary initial conditions. The six schemes are 
also compared by integrating them numerically from an initial state at 
rest with a height perturbation at a single point. The finite difference 
and finite element primitive equation schemes with unstaggered grid points 
give very poor results for the small scale features. The staggered scheme 
B gives much better results with both finite differences and finite 
elements. The vorticity-divergence system with unstaggered points also is 
very good with finite differences and finite elements. It is especially 
important to take into account these results when formulating efficient 
finite element prediction models. 
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1. IntFoduction 



The finite element method (FEM), which was developed in engineeping 
statics, has recently been introduced into various atmospheric prediction 
models (Cullen, 197^; Hinsman, 1975; Staniforth and Mitchell, 1977). The 
FEM is a special case of the Galerkin procedure in which the dependent 
variables are approximated by a finite sum of spatially varying basis 
functions with time dependent coefficients. The FEM basis functions are 
low order polynomials which are zero except in a localized region. The 
Galerkin procedure produces a set of coupled ordinary differential 
equations for the coefficients which are solved by introducing finite 
differences in time (see for example Finder and Gray (1977)). 

FEM models are potentially more accurate than finite difference 
models, but they normally require more computational effort per degree of 
freedom. For this reason it is especially important to formulate FEM 
models efficiently. Kelley and Williams (1976) found considerable small 
scale noise in an FEM model of the shallow water equations in a channel 
which had all variables carried at the same nodal points. Winninghoff 
(1968), Arakawa and Lamb (1977) and Schoenstadt (1980) have demonstrated 
the advantages of spatial staggering of dependent variables in finite 
difference models. Also Staniforth and Mitchell (1977, 1978) have 
obtained excellent results with a vorticity-divergence FEM formulation. 
This paper will compare these FEM formulations by considering the 
geostrophic adjustment process with the linearized shallow water 
equations in one dimension. 
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2. Basic Equations 

The linearized shallow-water equations with no mean flow can be 
written : 



3u p 4 . rr - n 
■ 5 ^ - £v + g ^ - 0 , 



3x 

+ fu = 0 , 
3u 



3v 

It 

^ + H = 0 , 
8t 9x 



( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 



where u and v are the perturbation velocities in the x and y directions, 
respectively, and H and h the mean and perturbed heights of the free 
surface. Also g represents gravity and f is the coriolis parameter. Note 
that all quantities are independent of y. 

The vorticity and divergence equations are obtained by differentiating 
(2.1) and (2.2) with respect to x which yields: 



3D , 3^h „ 



3C 

IT 



+ 



fD = 0 



« 



(2.4) 

(2.5) 



3t 



+ HD = 0 , 



( 2 . 6 ) 



where D = 9u/9x is the divergence and 4 = 9v/9x is the vorticity. These 
relations for D and ^ are particularly simple in this case since 9u/9y = 
9v/9y = 0. 

Schoenstadt (1977) solved the continuous equations (2.1)-(2.3) with 
the of the spatial Fourier transform. If we denote Fourier transforms by 
a tilde, such as 
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oo 



u(k,t) 



u(x,t) 



-ikx 

e 



dx 



..jOO 



(2.7) 



then the set (2.1)-(2.3) can be transformed to the form: 



^ = nfv - iygh , 
at 



( 2 . 8 ) 



dv 

dt 



= -nfu , 



(2.9) 



dh 

dt 



= -iyHu , 



( 2 . 10 ) 



where 0 = 1 and y = k. The quantities y and y will be useful later when 
finite difference and finite element solutions are needed. The initial 
conditions are written 



Uq = u(k,0) 



OO 




-jOO 



u(x,0) 



-ikx 

e 



dx , 



( 2 . 11 ) 



with similar definitions for v and h . Schoenstadt (1977) solved the set 

o o 

(2.7)-(2.9) subject to initial conditions by the eigenvalue-eigenvector 
approach which gives: 



u(k, t) 
v(k,t) 



V 

o o 

u cos vt + nf — sin vt sin vt , 

o V V 



nf _ 

V 'o 



2 u 2.2 

u sin vt + — P- + ^ .. cos vt} V. 



h(k,t) = - 



2 2 
V V 

+ {1 _ cos vt} , 

V 



u sin vt - {1 - cos vt}v, 

VO . .2 



V 

+{ *^ ^ ^ ■ cos vt} h^ , 

V V 



( 2 . 12 ) 



(2.13) 



(2.14) 



where : 



2 

V 



2 2 2 
n f + y gH . 



(2.15) 
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The transformed vorticity-divergence set (2.4)-(2.6) is written: 
dD ~ 2 ~ „ 

^ - fC - d gh = 0 , (2.16) 

(2.17) 



f + f6 - 0 . 



dh 

- + HD = 0 , 



(2,18) 



2 2 

where y = k . The solution to this set, which can be obtained directly 
or by using D = iku and C = ikv in (2 . 13)-(2 . 15) , is given by: 



where: 



D(k,t) 

C(k,t) 



D cos vt + sin vt + sin Vt , 

o V 



fD 



V 



sin vt + [ ^ f ^ + — cos vt] C 
2 

- [1 - cos vt] gh^ , 

V 



HD 



h(k,t) = - 



O . Hfr- ....I." 

Sin vt r 11 - cos Vt] C 

V ^2 



V = 



2 2 
f + y gH 



V V 



(2.19) 



( 2 . 20 ) 



( 2 . 21 ) 

( 2 . 22 ) 



3. Finite Difference and Finite Element Solutions 

Schoenstadt (1980) carried out a general analysis of the solutions to 
(2.1)-(2.3) which allowed for spatially centered finite differences or 
finite elements. We will use the same method to compare certain finite 
difference and finite element solutions to systems (2.1)-(2.3) and (2.4)- 
(2.6). The various finite difference and finite element forms correspond- 
ing to (2.1)-(2.3) or (2.4)-(2.6) are given in the Appendix. Following 
Schoenstadt (1980) the Fourier transformed versions of the various 
numerical schemes for the equations (2.1)-(2.3) can be written in the 
following form: 
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(3.1) 



a(k) — = fB(k) V - iga(k) h , 
a(k) ^ = -fB(k) u , 
a(k) = -iHa(k) u . 



(3.2) 

(3.3) 



The functions a(k), B(k) and o(k) are given in Table I for the various 
schemes considered. This set can be put in the same form as (2.8)-(2.10) 
by dividing by a and by setting: 

ri = S/a and y = a/a . (3.4) 

In this case the frequency is given by 

= (S^f^ + a^gH)/a^ . (3.5) 

The solutions to set (3.D-(3.3) are given by (2. 12)-(2. 14 ) with the use 
of (3.4) and (3.5). 

Table I. Coefficients in primitive equations for various 



numerical schemes. 






Scheme 


a 


6 


a 


differential 


1 


1 


k 


A 


1 


1 


sin (kAx)/Ax 


B 


1 


1 


sin (kAx/2)/( Ax/2) 


FEM A 


(2+cos(kAx) )/3 


(2+cos(kAx) )/3 


sin (kAx)/Ax 


FEM B 


(2+cos(kAx) )/3 


(2+cos(kAx) )/3 


(5sin(kAx/2)+sin(3kAx/2) )/4Ax 
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The numerical schemes for the vorticity-divergence system (2.^)-(2.6) 



lead to the following transformed equations: 

dD cr ^2 ~ ^ 

a 3 - - afC - o gh = 0 , 
at 



a 4r + fctD = 0 , 

dt 



dh 



(3.6) 

(3.7) 



a — + HaD = 0 , 
dt 



(3.8) 

where ct(k) and o(k) are given in Table II. This set can be put in the 
same form as ( 2 . 16 )- ( 2 . 18 ) by dividing by a and setting: 



2 ^2 , 
y = o /a 



(3.9) 



The frequency equation (2.22) becomes 



+ (o^/a)gH 



(3.10) 



which has a different form from (3.5). The solutions to set (3.6)-(3.8) 
are given by ( 3 . 6 )-( 3 . 8 ) with the use of ( 3 . 9 ) and ( 3 . 10 ). 



Table II. 

Scheme 

Scheme 

differential 



Coefficients 
equations for 
a 



in vorticity-divergence 
various numerical schemes. 



1 



k 



2 



finite difference 



sin2(kAx/2)/(Ax/2)^ 



FEM 



(2+cos(kAx ) ) /3 sin2(kAx/2)/( Ax/2) 2 
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The various parameters which determine the solutions (2. 13)-(2. 15 ) and 
(2. 19)-(2.21) are shown in Tables I and II, respectively. Table I contains 
Schemes A and B for the primitive equations where Scheme A is unstaggered 
and Scheme B has the velocity points midway between the height points (see 
Schoenstadt, 1980). The table also includes the finite element forms 
which are obtained when piecewise linear basis functions are used. Note 
that k is poorly represented by o with Scheme A near k = tt/Ax, and that 
the problem remains with the FEM version of Scheme A. The staggered grid 
gives a much better approximation since spatial derivatives are computed 
over a distance of Ax compared to 2Ax with the unstaggered grid. 

Table II contains the parameters for the finite difference and finite 
element versions of the vorticity-divergence set of equations. In this 

case vorticity, divergence and height are carried at the same points. 

2 2 
Note that o for both cases is the same as the value of a for Scheme B 

from Table I. It can be seen from the tables that the staggered primitive 

equation (Scheme B) and vorticity-divergence formulations have the same 

values for a and a and therefore for v, so that these should give the same 

solution except for truncation error in the initial conditions. 

As pointed out by Schoenstadt (1980), the solutions (2. 12)-(2. 1^1 ) for 
the various schemes differ only through the coefficients H/v, PA*, and 
nP/v , and the same dependence occurs in system (2. 19 )-(2. 21 ) with n = 1, 
except that the coefficient np/v^ does not appear. Figure la contains 
the phase velocity, c = v/k , as a function of kAx/rr for the various 
schemes in Tables I and II as computed from (3.5) and (3.10), respec- 
tively. The differential solution approaches f/k for small k and 
the shallow-water speed (gH)^^^ for large k. Scheme A gives the poorest 
phase speed and the finite element Scheme A is also very poor for the 
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highest wavenumbers. The finite element scheme B is very close to the 

differential solution, while the vorticity-divergence FEM scheme is a 

little higher. The group velocity, G = dv/dp, is given in Fig. lb, as a 

function of kAx/ir. The differential solution is zero at k = 0 and it 

1 /2 

approaches the shallow-water phase speed (gH) for large k. Scheme A 
and its FEM version are very poor for the short waves since they 
propagate energy in the wrong direction. In fact the FEM scheme gives a 
group velocity which is more than double the correct value and of the 
wrong sign, at certain points. The FEM scheme B gives the best group 
velocity while the FEM vorticity-divergence scheme gives values that are 
somewhat higher. 

2 

The coefficients n/v, y/v and ny/v are given in Fig. 2a, 2b and 2c, 
respectively, as functions of kAx/7T. Scheme A is the poorest for each 
coefficient, but the FEM version of scheme A is just as bad for the short 
waves. The best scheme is the FEM version of scheme B, although the FEM 
vorticity-divergence scheme is also very good. The coefficient n/v, which 
is given in Fig. 2a, is especially important since n /v relates the 
initial height to the final (steady-state) height field (see (2.14)). In 
particular, the figure shows that if = 0, the final h for k = tt/Ax is 

more than 25 times too large for scheme A and the FEM version of A! This 
is one reason why non-staggered schemes tend to generate small scale 
noise. These results were given by Schoenstadt (1980) with the exception 
of the vorticity-divergence schemes. 

4. Final State Example 

The two aspects of the geostrophic adjustment process that must be 
considered in assessing a particular numerical scheme are: 1) forecast 
time required to reach the adjusted state, 2) the accuracy of the final 
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adjusted state. The group velocity curves in Fig. la provide an 
indication of the comparative adjustment times for the various schemes. 

The final adjusted state, which is more important, could be obtained by 
Fourier transforming the terms that are independent of t in (2. 12)-(2. 14) 
or (2. 19 )-(2. 21 ) . However, in this paper the final state will be 
determined by integrating the finite difference equations in t until the 
adjusted state is reached. This approach is preferable because time 
differencing effects are included and a time filter can also be used. 

The various sets of equations, which are given in the Appendix are 
integrated with centered time differences. The time filter developed by 
Robert (1966) (see also Asselin, 1972) is applied to the past time value 
with the coefficient y = .05. The new time values for the FEM schemes 
are found by Gauss elimination. 

The initial conditions are given by: 

a I X I £ Ax/ 2 

h(x,o) = { (4.1) 

o |x| > Ax/2 , 

u(x,o) = v(x,o) = 0 , or C(x,o) = D(x,o) = 0 . 

These initial conditions are convenient for comparing the various schemes 
since no truncation error is introduced when the initial vorticity and 
divergence are computed from these initial velocities. The analytic 
solution for the final adjusted h field can be obtained by integrating the 
following expression that was obtained by Schoenstadt (1977): 

oo 

h (x) = h(x,o) ^ sgn(x-C)e~l^~^l^^[-| |^(^,0)-v(?,0)]dC , 

..JOO 

(4.2) 
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1 /2 

where h^(x) is the final adjusted height and X = (gH) /f is the Rossby 
radius of deformation. The initial geostrophic wind which is required in 
(4.2) can be conveniently written: 

I ||(x,0) = ^ [6(x+Ax/2)-6(x-Ax/2)] , (4.3) 

where 6(x) is the delta function. 

When (4.1) and (4.3) are introduced into (4.2) the solution becomes: 

e sinh(Ax/2A) Ax/2 < x 

h (x) = a 1 - e cosh(x/A) -Ax/2 < x < Ax/2 . (4.4) 

s — — 

e^^^ sinh(Ax/2A) x < -Ax/2 

Fig. 3 contains h^(x) for the case Ax = A/2. 

The numerical integrations with the various schemes are performed on a 
grid of 200 points with cyclic boundary conditions. The initial distur- 
bance at X = 0 is placed in the center of the computational domain so that 
the cyclic boundary conditions will not affect the solution near x = 0 
until well after the adjusted state is reached. Fig. 3 includes the numer- 
ical solutions at t = 3 days for the following schemes: A, B and FEM A. 

Scheme A shows strong oscillations with every other point returning to 0. 
The FEM scheme A has smaller oscillations near x = 0, but they become 
larger than the oscillations with scheme A farther out. Scheme B gives 
very smooth behavior and it is close to the analytic solution. The 
vorticity-divergence system gives the same solution as scheme B, and is 
very close to the analytic solution as can be seen in Table III which 
compares the solutions at the first two grid points. 
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Table III. Numerical solutions h/a at t = 72 hours for the first two grid 
points for various schemes compared with analytic solution. 



X 


0 


Ax 


Differential 


0.221 


0. 153 


A 


0.459 


0.0 


B 


0.240 


0.148 


vorticity-divergence 


0.240 


0. 148 


FEM A 


0.298 


0.084 


FEM B 


0.227 


0. 157 


FEM vorticity-divergence 


0.213 


0. 154 



The results given in Fig. 3 and Table III are consistent with the 

2 2 

curves for h /v shown in Fig. 2a, since h is proportional to n /v (see 

s 

(2.14) and (2.21)). In particular the poor behavior for the unstaggered 
primitive equation schemes (A and FEM A) in Fig. 2a is consistent with the 
large amplitude short waves in Fig. 3. Also the large oscillations 
farther out with FEM A may be the result of the large spurious group 
velocity that is shown in Fig. 1b for that scheme. All the staggered 
primitive equation and vorticity-divergence schemes give excellent 
predictions of the final adjusted height field. It should be pointed out 
that the inclusion of light time smoothing (Y = .05) is necessary to 
produce the spatially smooth solutions for these cases. Apparently the 
vanishing group velocity for kAx/ir = 1 (see Fig. 1b) does not allow the 
smallest scale gravity waves to propagate out from the initial distur- 
bance. Haltiner and McCollough (1975) demonstrated the usefulness of time 
filtering in a baroclinic primitive equation model. 
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5. Conclusions 



The objective of this paper is to determine the response of various 
finite element schemes to small scale initial conditions or small scale 
forcing. It is especially important that FEM prediction schemes properly 
describe small scale features, because FEM models usually require more 
computational effort per degree of freedom than most finite difference 
models. This study treated the geostrophic adjustment process with the 
linearized primitive equations and also with the related vorticity- 
divergence set of equations. The development followed Schoenstadt (1980) 
wherein the spatially discretized equations were Fourier transformed in x, 
and then solved with arbitrary initial conditions. These solutions were 
dependent on certain coefficients which were computed for the various 
numerical schemes and compared with the differential expressions. Three 
FEM schemes were examined as well as the three corresponding finite 
difference schemes. It was found that the unstaggered (scheme A) primi- 
tive equation model gives the poorest behavior followed by the correspond- 
ing FEM formulation. These schemes are especially bad for the shortest 
resolvable scales. The finite difference primitive equation model, which 
staggers height points between velocity points (scheme B) has much better 
behavior than the unstaggered schemes. The vorticity-divergence model 
where C, D and h are carried at the same points has the same coefficients 
as scheme B. The FEM version of scheme B, which has staggered nodal points, 
was found to have the best behavior and the FEM vorticity-divergence model 
was also found to be very good. 

The six schemes were also compared by integrating them numerically with 
centered time differences from an initial state at rest with a height 
perturbation at a single point. The analytic solution for this initial 
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state approached a smooth height field after the inertial gravity waves 
radiated away. Scheme A and the FEM form of scheme A gave very poor 
solutions with large oscillations from point to point. All of the other 
schemes produced smooth solutions with the FEM schemes being the most 
accurate. The smoothness of these solutions was improved by light time 
smoothing. Although the initial state used in this comparison is somewhat 
extreme, it shows clearly the superiority of the staggered primitive 
equation and vorticity-divergence schemes over the non-staggered primitive 
equation schemes. 

Winninghoff (1968), Arakawa and Lamb (1977) and Schoenstadt (1980) 
have demonstrated the advantages of spatial staggering of predictive 
variables in finite difference models. Our results strongly indicate that 
FEM models should either use staggered nodal points in the primitive equations 
or unstaggered nodal points in the vorticity-divergence equations (see 
also Schoenstadt, 1980). In fact Staniforth and Mitchell (1977, 1978) 
have developed a FEM model based on the vorticity-divergence form of the 
shallow-water equations that produces smooth forecasts with only time 
smoothing. In contrast, Kelley and Williams (1976) obtained very noisy 
results with an unstaggered FEM model which used the primitive equations 
for flow in a channel. If non-staggered finite FEM element models are 
used, it is often necessary to use high order smoothing to damp the small 
scales as discussed by Cullen (1976). Thacker (1978) tested a finite 
element formulation of the linearized shallow-water equations with 
staggered nodal points and he obtained smooth solutions. 

Since FEM models usually require more computer time per degree of 
freedom, it is very important for the numerical scheme used to be accurate 
for as small a scale as possible. In this paper we have shown that the 
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usual non-staggered FEM formulation of the primitive equations gives very 
poor geostrophic adjustment for small scale initial conditions. The same 
conclusion follows for small scale heating. On the other hand either the 
use of the primitive equations with staggered nodal points or the 
vorticity-divergence equations with unstaggered nodal points gives 
excellent treatment of small scale features in the geostrophic adjustment 
process. Clearly, the use of either formulation should be much more 
efficient than the unstaggered primitive equations, even when the latter 
have smoothing to destroy the smallest scale features. 
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Appendix 

In this Appendix the spatially discretized prediction equations are 
given for each scheme with x = mAx. The following schemes approximate the 
primitive equations (2.1)-(2.3): 
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FEM Scheme B 
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where Mot = (a . + 4a +a ,)/6. 
m m + 1 m m- 1 

Scheme B is staggered in such a way that the height points are 
equi-distant between the velocity points. The FEM equations can be 
derived with piecewise linear basis functions (see for example Chapter 6 
in Haltiner and Williams, 1980). 

The vorticity-divergence system (2.4)-(2.6) is approximated with the 
following schemes. 

Vorticity-Divergence 
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Finite Element Vorticity-Divergence 
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Figure Captions 



Fig. 1 



Fig. 2 



Fig. 3 



. The phase velocity c = v/y, and the group velocity G = dv/dy as 
functions of IcAx/tt for various numerical schemes. The curves 
are labeled as follows: 1) differential solution, 2) scheme A, 

3) scheme B and vorticity divergence finite difference scheme, 

4) FEM scheme A, 5) FEM scheme B, 6) FEM vorticity -divergence 
scheme.. These results use the following values: 

gH = 10^ f = Ax = 500 km. 

2 

. The coefficients r|/v, y/v and ny/v as functions cf kAx/7T, with 
the same labeling as in Fig. 1. 



The numerical solutions for schemes A, B and FEM. A as functions 
of x/Ax at t = 3 days. The steady-state differential solution, 
which is given by (4.4), is included for comparison. 
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